close all
clear
clc
%%%%%%%%%%%%%%%%%%%%%
%     Parameter     %
%%%%%%%%%%%%%%%%%%%%%
a=0.2;                      %lowpass filter not used
slideWindow=40;             %delta function forwarding and tailing window
cutoff=1.8;                  %peak selection cutoff critias
excludePercentage=0;      %fitting excluding percentage
fitWindow=10;
cutoffForce=12.1;
stepCutoffMagnet=12.1;
internalTTest=1;            %Single T test
internalTTestAlpha=0.01;
position_overstretching=12.41;
aa_length=0.5;    % the pesistent length of a amino acid
filelist=importdata('G:\Mingxi\Dropbox\yingjian guo\Alphacat\filelist-jump.txt');
only_unfolding=1;% only process the data in the unfolding part of spectrum
decomposecycles=0; % Separate data into different forcecycles by addition of a numerical flag (only the force altering region is outpeed
only_refolding=0;
with_smoothing=0;%whether to smooth rawdata
drift_correction=1; %whether to correct linear drift in long time measurements
forcechange_correction=0;%whether force is changing when steps are detected changes fitting parameters.
forcechange_correction_direct=1;
%change to dump if the analysis is repeated to avoid double counting
summaryoutputfile=('G:\Mingxi\Dropbox\yingjian guo\Alphacat\goodunfolding.txt');
good_threshold=1;


tic
index=1;
for currentfile = filelist.textdata'
    
    name=char(currentfile);
    
    
    %% Data processing
    steps_data=[];
    cycledata=[];
    filename=sprintf('%s%s',name,'.txt')
    %c=65/(exp(-(13.5-position_overstretching)/0.36)+0.48*exp(-(13.5-position_overstretching)/1.12));
    c=filelist.data(index)
    index=index+1;
    %     datastruct=importdata(filename);
    %     dataraw=datastruct.data;
    %     egxT=dataraw(:,1)';
    %     egMagneticT=dataraw(:,2);
    %     egRawDynT=dataraw(:,3);
    
    
    % Open file
    fid3 = fopen(filename); % Open file
    
    dataraw=textscan(fid3,'%f %f %f %f %f %f %f %f %f %f','HeaderLines',2);
    % Close file
    fclose(fid3);
    
    
    %     datastruct=importdata(filename);
    %     dataraw=datastruct.data;
    egxT=cell2mat(dataraw(1))';
    egMagneticT=cell2mat(dataraw(2));
    egRawDynT=cell2mat(dataraw(3));
    
    %plot(egxT); % Check equal spacing
    
    %Mingxi data: A. linear drift remover
    %             B. force altering regions
    %             C. points for fitting selection
    
    if drift_correction
        linearFitP=polyfit(egxT(egMagneticT==min(egMagneticT))',egRawDynT(egMagneticT==min(egMagneticT)),1);
        egfittedDriftT=polyval(linearFitP, egxT)';
        disp('Determining the drift\n');
        linearFitP % display the drift
        egRawDynTShift=egRawDynT-egfittedDriftT;
    else
        egRawDynTShift=egRawDynT;
    end
    %figure;
    %plot(egxT,egRawDynT);
    %hold on; plot(egxT,egfittedDriftT,'r','LineWidth',3);
    %figure;
    %plot(egxT,egRawDynT-egfittedDriftT,'b');
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %cutoffForce=min(egMagneticT);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %process if only unfolding signal is processed
    tsig=egMagneticT>cutoffForce;
    dsig=diff([0;tsig;0]);
    startIndex=find(dsig > 0);
    endIndex=find(dsig < 0)-2;
    
    
    %     hFig = figure();
    %     set(hFig, 'Visible', 'off')
    
    %uses parallel for to increase speed,need some testing, to do it serial
    %just change back to for, disabled ploting due to incompability to
    %parfor
    
    parfor forceAlteringRegionIndex=1:length(startIndex)
        is_good=0;% whether this tether has enough unfolding steps
        single_steps_data=[];
        single_cumulative_contour=0; % total unfolding size in a single force altering region
        region_end_index=endIndex(forceAlteringRegionIndex);
        region_start_index=startIndex(forceAlteringRegionIndex);
        if only_unfolding    %with this option only the signals in the unfolding region of the trajectory is analyzed
            region_Magnetic=egMagneticT(startIndex(forceAlteringRegionIndex):endIndex(forceAlteringRegionIndex)); % find the magnetic positions in the region
            region_magsig=diff([0;region_Magnetic;0]);
            if ~isempty(find(region_magsig<0,1,'first'))
                region_end_index=region_start_index+find(region_magsig<0,1,'first')-1;
            end
        end
        
        if only_refolding    %with this option only the signals in the refolding region of the trajectory is analyzed
            region_Magnetic=egMagneticT(startIndex(forceAlteringRegionIndex):endIndex(forceAlteringRegionIndex)); % find the magnetic positions in the region
            region_magsig=diff([0;region_Magnetic;0]);
            if ~isempty(find(region_magsig<0,1,'first')) && find(region_magsig<0,1,'first')<region_end_index
                region_start_index=region_start_index+find(region_magsig<0,1,'first');
            end
        end
        
        %find the postion where unfolding and refolding shifts
        if decomposecycles && ~only_refolding && ~only_unfolding
            region_Magnetic=egMagneticT(startIndex(forceAlteringRegionIndex):endIndex(forceAlteringRegionIndex)); % find the magnetic positions in the region
            region_magsig=diff([0;region_Magnetic;0]);
            if ~isempty(find(region_magsig<0,1,'first')) && find(region_magsig<0,1,'first')<region_end_index
                force_revert_position=find(region_magsig<0,1,'first')-1;
            end
        end
        
        %region_start_index
        %region_end_index
        egx=egxT(region_start_index:region_end_index);
        egRawDyn=egRawDynTShift(region_start_index:region_end_index);
        %scaled magnetic Location
        egMagnetic=egMagneticT(region_start_index:region_end_index);
        
        %generate flagto mark the unfolding and refolding cycles
        if decomposecycles
            if  only_unfolding
                egcycleflag=forceAlteringRegionIndex.*ones(size(egx));
            elseif only_refolding
                egcycleflag=(-forceAlteringRegionIndex).*ones(size(egx));
            else
                egcycleflag=(-forceAlteringRegionIndex).*ones(size(egx));
                egcycleflag(1:force_revert_position)=forceAlteringRegionIndex;
            end
        end
        %egTruDyn=(egMagnetic-min(egMagnetic))*(max(egRawDyn)-min(egRawDyn))/(max(egMagnetic)-min(egMagnetic));
        %figure;
        %plot(egxT,egRawDynTShift);
        %hold on;plot(egx,egRawDyn,'r');
        
        
        %% Assign
        x=egx;
        RawDyn=egRawDyn;
        
        
        %TruDyn=egTruDyn;
        %% Filtering By lowpass filter (Freqency**)
        % Without filter seems better for loading rate data (or change parameter)
        FilDynA=filter(a,[1 a-1],RawDyn);
        FilDynB=smooth(RawDyn,slideWindow,'sgolay'); % SG smooth
        %
        
        %         subplot(3,1,1);
        %         plot(x,RawDyn,'b',x,FilDynA,'r');
        %         hold on; plot(x,FilDynB,'g');
        
        
        %% auxiliary function (delta function)
        
        disp([num2str(slideWindow), ' points delta function ...']);
        if with_smoothing
            customDyn=FilDynB;
        else
            customDyn=RawDyn;
        end
        
        %Apply high pass filter to get rid of slow increases of force by a
        %linear function
        
        if forcechange_correction
            ignore_points=100;
            loadingjumpfit=polyfit((1:(size(RawDyn)-ignore_points))',RawDyn((ignore_points+1):end),1);
            deltaDyn=RawDyn -(loadingjumpfit(1).*(1:(size(RawDyn)))');
        end
        
        if forcechange_correction_direct
            ignore_points=100;
            loadingjumpfit=polyfit((1:(size(RawDyn)-ignore_points))',RawDyn((ignore_points+1):end),1);
            deltaDyn=RawDyn -(loadingjumpfit(1).*(1:(size(RawDyn)))');
            RawDyn=deltaDyn;
        end
        %increment the cycle outputdata
        if decomposecycles
            egforce=c.*(exp(-(13.5-egMagnetic)./0.36)+0.48.*exp(-(13.5-egMagnetic)./1.12));
            cycledata=[cycledata;[x',egMagnetic,egforce,customDyn,egcycleflag']];
        end
        deltaF=zeros(1,numel(customDyn));
        
        %starting
        if forcechange_correction
            for i=1:slideWindow
                diff_value=i-1;
                forwarding=mean(deltaDyn(i:i+diff_value));
                tailing=mean(deltaDyn(i-diff_value:i));
                deltaF(i)=tailing-forwarding;
            end
            
            for i=slideWindow+1:numel(deltaDyn)-slideWindow
                %disp(i);
                forwarding=mean(deltaDyn(i:i+slideWindow));
                tailing=mean(deltaDyn(i-slideWindow:i));
                deltaF(i)=tailing-forwarding;
            end
            
            %ending
            for i=numel(deltaDyn)-slideWindow+1:numel(deltaDyn)
                diff_value=numel(deltaDyn)-i;
                forwarding=mean(deltaDyn(i:i+diff_value));
                tailing=mean(deltaDyn(i-diff_value:i));
                deltaF(i)=tailing-forwarding;
            end
        else
            for i=1:slideWindow
                diff_value=i-1;
                forwarding=mean(customDyn(i:i+diff_value));
                tailing=mean(customDyn(i-diff_value:i));
                deltaF(i)=tailing-forwarding;
            end
            
            for i=slideWindow+1:numel(customDyn)-slideWindow
                %disp(i);
                forwarding=mean(customDyn(i:i+slideWindow));
                tailing=mean(customDyn(i-slideWindow:i));
                deltaF(i)=tailing-forwarding;
            end
            
            %ending
            for i=numel(customDyn)-slideWindow+1:numel(customDyn)
                diff_value=numel(customDyn)-i;
                forwarding=mean(customDyn(i:i+diff_value));
                tailing=mean(customDyn(i-diff_value:i));
                deltaF(i)=tailing-forwarding;
            end
        end
        
        %% peak location determination
        
        disp([num2str(cutoff), ' cutoff peak locating ...']);
        
        deltaFCut=deltaF.*(deltaF>=cutoff|deltaF<=-cutoff);
        deltaFCutU=deltaF.*(deltaF>=cutoff);  %positive
        deltaFCutD=deltaF.*(deltaF<=-cutoff); %negative
        
        
        [Upeaks,Ulocs]=findpeaks(deltaFCutU,'minpeakdistance',slideWindow);
        [Dpeaks,Dlocs]=findpeaks(-deltaFCutD,'minpeakdistance',slideWindow);
        
        %         subplot(3,1,3);
        %         plot(x,deltaF,'k-');
        %         hold on;plot(x,deltaFCut,'g');
        %         plot(Ulocs*(max(x)-min(x))/(numel(x)-1)+min(x),Upeaks,'ro');
        %         plot(Dlocs*(max(x)-min(x))/(numel(x)-1)+min(x),-Dpeaks,'mo');
        
        %%%%%%%%%%%%%%%%%%%%%%
        %        TRAIL       %
        %%%%%%%%%%%%%%%%%%%%%%
        %% T Test filter remove false steps
        fitDyn=RawDyn;
        locs=sort([1,Ulocs, Dlocs, numel(fitDyn)]);
        
        % Multiple rounds of test until fix
        if forcechange_correction
            if internalTTest == 1
                compareDyn=deltaDyn;
                compareDyn_size=size(compareDyn);
                templocs=zeros(1,numel(locs));
                previouslocs=0;
                count=1;
                
                while numel(previouslocs) ~= numel(locs)
                    disp([num2str(internalTTestAlpha*100), '% Internal Student T Test: ', num2str(count)]);
                    previouslocs = locs;
                    for i=2:numel(locs)-1
                        %interval selection
                        startI=locs(i)-fitWindow;
                        endI=locs(i)+fitWindow;
                        if startI < 1
                            startI=1;
                        end
                        if endI > compareDyn_size
                            endI=compareDyn_size;
                        end
                        
                        jumpI=locs(i);
                        startCutI=startI+floor((jumpI-startI)*excludePercentage);
                        endCutI=endI-floor((endI-jumpI)*excludePercentage);
                        groupA=compareDyn(startCutI:jumpI-1);
                        groupB=compareDyn(jumpI+1:endCutI);
                        
                        [H,PI]=ttest2(groupA,groupB,internalTTestAlpha,[],'unequal');
                        %disp(PI);
                        
                        if H==1
                            templocs(i)=locs(i);
                        else
                            locs(i)=startI;
                        end
                    end
                    
                    locs=[1,templocs(templocs~=0), numel(fitDyn)];
                    templocs=zeros(1,numel(locs));
                    count=count+1;
                end
            end
        else
            if internalTTest == 1
                compareDyn=fitDyn;
                compareDyn_size=size(compareDyn);
                templocs=zeros(1,numel(locs));
                previouslocs=0;
                count=1;
                
                while numel(previouslocs) ~= numel(locs)
                    disp([num2str(internalTTestAlpha*100), '% Internal Student T Test: ', num2str(count)]);
                    previouslocs = locs;
                    for i=2:numel(locs)-1
                        %interval selection
                        startI=locs(i)-fitWindow;
                        endI=locs(i)+fitWindow;
                        if startI < 1
                            startI=1;
                        end
                        if endI > compareDyn_size
                            endI=compareDyn_size;
                        end
                        
                        jumpI=locs(i);
                        startCutI=startI+floor((jumpI-startI)*excludePercentage);
                        endCutI=endI-floor((endI-jumpI)*excludePercentage);
                        groupA=compareDyn(startCutI:jumpI-1);
                        groupB=compareDyn(jumpI+1:endCutI);
                        
                        [H,PI]=ttest2(groupA,groupB,internalTTestAlpha,[],'unequal');
                        %disp(PI);
                        
                        if H==1
                            templocs(i)=locs(i);
                        else
                            locs(i)=startI;
                        end
                    end
                    
                    locs=[1,templocs(templocs~=0), numel(fitDyn)];
                    templocs=zeros(1,numel(locs));
                    count=count+1;
                end
            end
        end
        
        %% least square fit to heaviside function (per interval)
        %Change to fit locally to adjacentpoitns
        fittedDyn=zeros(numel(fitDyn),1);
        
        for i=2:numel(locs)-1
            %tempfix for some locs<0 error
            if locs(i)>fitWindow && locs(i)<(length(x)-fitWindow)
                %interval selection
                startI=locs(i)-fitWindow;
                endI=locs(i)+fitWindow;
                jumpI=locs(i);
                startCutI=startI+floor((jumpI-startI)*excludePercentage);
                endCutI=endI-floor((endI-jumpI)*excludePercentage);
                currentX=x(startCutI:endCutI);
                currentX_index=startCutI:endCutI;
                currentY=fitDyn(startCutI:endCutI);
                
                %figure;plot(x(startI:endI),fitDyn(startI:endI),'go',currentX,currentY,'ro')
                %guess starting points
                ga=mean(fitDyn(jumpI:endCutI))-mean(fitDyn(startCutI:jumpI));
                gb=jumpI;
                if ga>=0
                    gc=mean(fitDyn(startCutI:jumpI));
                else
                    gc=mean(fitDyn(jumpI:endCutI));
                end
                
                %fit (added a linear term to approximate the loading rate
                if forcechange_correction
                    gd=polyfit((startCutI:jumpI)',fitDyn(startCutI:jumpI),1);
                    d=gd(1);
                    fequation='a.*heaviside(x-b)+c+d.*(x-b)';
                    ft = fittype(fequation,'problem',{'d'});
                    fo = fitoptions('Method','NonlinearLeastSquares','StartPoint',[ga gb gc],'Lower',[-inf 0 -inf],'Upper',[inf inf inf]);
                    currentFit = fit(currentX_index',currentY,ft,fo,'problem',d);
                    jumptime=x(locs(i))-min(x);
                    locs(i)=round(currentFit.b);
                    %                 leftValue=currentFit.c;
                    %                 rightValue=currentFit.a+currentFit.c;
                    %                 stepValue=0.5*currentFit.a+currentFit.c;
                    
                else
                    fequation='a.*heaviside(x-b)+c';
                    ft = fittype(fequation);
                    fo = fitoptions('Method','NonlinearLeastSquares','StartPoint',[ga gb gc],'Lower',[-inf 0 -inf],'Upper',[inf inf inf]);
                    currentFit = fit(currentX_index',currentY,ft,fo);
                    jumptime=x(locs(i))-min(x);
                    locs(i)=round(currentFit.b);
                    %                 leftValue=currentFit.c;
                    %                 rightValue=currentFit.a+currentFit.c;
                    %                 stepValue=0.5*currentFit.a+currentFit.c;
                end
                
                
                %currentF=currentFit(currentX);
                %hold on;plot(currentX,currentF,'k','LineWidth',3);
                %currentFit
                
                %update result by directly fitting index
                %locs(i)=round((currentFit.b-min(x))/(max(x)-min(x))*(numel(x)-1));
                
                %is this the time from starting point?
                
                %                 if locs(i)>fitWindow
                %                     fittedDyn((locs(i)-fitWindow):locs(i)-1)=leftValue;
                %                 else
                %                     fittedDyn(1:locs(i)-1)=leftValue;
                %                 end
                %                 fittedDyn((locs(i)+1):locs(i)+fitWindow)=rightValue;
                %                 fittedDyn(locs(i))=stepValue;
                %get the parameters for the magnet and calculated the force
                stepsize=currentFit.a;
                locmagnet=egMagnetic(locs(i));
                force_then=c*(exp(-(13.5-locmagnet)/0.36)+0.48*exp(-(13.5-locmagnet)/1.12));
                full_length=stepsize/zz_wlc(force_then*aa_length/4);
                %record data for forces > preset values to reduce false positives
                if locmagnet>stepCutoffMagnet
                    steps_data=[steps_data;[forceAlteringRegionIndex,locmagnet,force_then,stepsize,jumptime]];
                    single_cumulative_contour=single_cumulative_contour+full_length;
                    single_steps_data=[single_steps_data;[forceAlteringRegionIndex,locmagnet,force_then,stepsize,jumptime]];
                end
            end
            
        end
        
        
        if single_cumulative_contour>=good_threshold
            is_good=1;
        end
        
        %write to global file if the tether is good
        if is_good
            fid1=fopen(summaryoutputfile,'a');
            for row= 1:size(single_steps_data,1)
                fprintf(fid1,'%s\t%d\t%.3f\t%.3f\t%.3f\t%.3f\n',name,single_steps_data(row,:)');
            end
            fclose(fid1);
        end
        
        % %update last step & first value
        % rightValue=currentFit.c+currentFit.a;
        % fittedDyn(locs(i)+1:numel(fittedDyn))=rightValue;
        % fittedDyn(1)=fittedDyn(2);
        
        % %update locs value
        % for i=2:numel(locs)-1
        %     fittedDyn(locs(i))=.5*(fittedDyn(locs(i)-1)+fittedDyn(locs(i)+1));
        % end
        
        %         subplot(3,1,2);
        %         plot(x,fitDyn,'ro');
        %         hold on; plot(x,fittedDyn,'g','LineWidth',1);
        %         figure_filename = sprintf('%s%s',name,'.eps');
        %         print(gcf, '-depsc2','-painters', figure_filename);
        %
    end
    
    
    %open write file
    if only_unfolding
        filename = sprintf('%s%s',name,'-f-unfoldingstep-jump.dat');
    elseif only_refolding
        filename = sprintf('%s%s',name,'-f-refoldingstep-jump.dat');
    else
        filename = sprintf('%s%s',name,'-f-step.dat');
    end
    fid = fopen(filename,'w');
    if(fid == 0)
        disp('file cannot be opened');
    end
    fprintf(fid,'%s\n',name);
    fprintf(fid,'rampnum\tmagnet(mm)\tforce (pN)\tstepsize(nm)\tconter length(nm)\n');
    fprintf(fid,'%d\t%.3f\t%.3f\t%.3f\t%.3f\n',steps_data');
    fclose(fid);
    
    %open the file for input the data and marking for the cycles
    
    if decomposecycles
        if only_unfolding
            cyclefilename = sprintf('%s%s',name,'-cycle-unfoldingdata.dat');
        elseif only_refolding
            cyclefilename = sprintf('%s%s',name,'-cycle-refoldingdata.dat');
        else
            cyclefilename = sprintf('%s%s',name,'-cycle-data.dat');
        end
        
        fid2 = fopen(cyclefilename,'w');
        if(fid2 == 0)
            disp('file cannot be opened');
        end
        fprintf(fid2,'%s\n',name);
        fprintf(fid2,'Time\tmagnet(mm)\tforce (pN)\tExtension(nm)\tcycleflag(nm)\n');
        fprintf(fid2,'%d\t%.3f\t%.3f\t%.3f\t%.3f\n',cycledata');
        fclose(fid2);
    end
    %open force cycle decompos
    
    
end
toc